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ABSTRACT 


A  procedure  for  investigating  glass  cladding  behavior  under  arbitrary  loads, 
including  fluctuating  wind  loads,  is  presented.     The  procedure  accounts  for 
the  fact  that  internal  stresses  are  nonlinear  functions  of  the  external  loads 
that  initial  glass  strengths  are  random  functions  of  position  and  direction, 
and  that  glass  strength  undergoes  degradation  under  the  action  of  external 
loads  in  accordance  with  basic  fracture  mechanics  laws.     Numerical  examples 
are  presented,  and  corresponding  probability  distribution  curves  are  calculat 
indicating  the  probability  of  failure  of  a  specified  panel  subjected  to  fluc- 
tuating wind  loads  and  to  1-minute  constant  loads.     These  curves  are  used  to 
illustrate  a  method  for  assessing  current  glass  cladding  design  procedures. 
For  the  case  considered  in  the  paper,  it  was  found  that  transformation  of 
the  peak  wind  load  averaged  over  1-2  seconds  into  an  equivalent  1-minute  load 
appears  to  underestimate  the  probability  of  failure  of  glass  cladding.  The 
work  reported  in  the  paper  is  part  of  an  ongoing  window  cladding  research 
program  being  conducted  at  the  National  Bureau  of  Standards. 

Key  words:     aerodynamics;  buildings;  deformation;  engineering  mechanics; 
failure;  glass;  loads  (forces);  probability  theory. 
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1.  INTRODUCTION 


The  improvement  of  design  procedures  for  cladding  glass  subjected  to  wind 
loads  has  been  the  focus  of  a  considerable  amount  of  research  in  recent  years. 
In  1979,  Pittsburgh  Plate  Glass  (PPG)  published  revised  design  charts  [2]  based 
upon  nonlinear  stress  analyses  used  in  conjunction  with  elementary  statistical 
methods.     In  a  number  of  instances,  these  charts  differ  from  those  issued  by 
other  manufacturers,  notably  LOF  [8],  and  exhibit  internal  inconsistencies  as 
well.     For  example,  according  to  the  PPG  charts,  a  3.67  ft  x  18.33  ft  (1.12  m 
X  5.60  m)  glass  panel  with  a  thickness  of  12.7  mm  supported  on  four  sides 
deflects  11.9  mm  under  a  load  of  100  psf  (4790  Pa));  for  a  panel  with  a  span 
of  1.12  m,  with  the  same  thickness  and  under  the  same  load  but  supported  on 
two  sides,  the  deflection  as  obtained  from  the  PPG  charts  is  only  9.1  mm, 
rather  than  being  equal  to  or  in  excess  of  11.9  mm.     Inconsistencies  with 
respect  to  design  loads  exist  as  well.     As  an  example,  for  a  1/8  in  (3  mm) 
thick  annealed  float  glass  panel  supported  on  four  sides  the  PPG  charts  specify 
a  15  psf  (719  Pa)  one-minute  load  if  the  dimensions  of  the  panel  are  2.83  ft  x 
7.07  ft  (0.865  m  X  2.15  m) ,  and  the  same  -  rather  than  a  larger  -  one-minute 
load  if  the  dimensions  are  2.24  ft  x  6.70  ft  (0.68  m  x  2.04  m) .  Inconsisten- 
cies such  as  these,  and  discrepancies  with  respect  to  other  manufacturers' 
charts  [8],  suggest  that  the  development  of  an  improved  theoretical  framework 
for  the  design  of  glass  cladding  is  a  necessary  task. 

An  important  step  toward  such  an  improvement  was  proposed  in  1980  by  Beason 
[2].     Reference  2  combined  nonlinear  stress  analysis  with  the  classic  Weibull 
theory  to  estimate  the  probability  of  failure  of  a  glass  panel  subjected  to  a 
specified  load,  given  the  parameters  of  the  Weibull  distribution  of  the  glass 
strength.     Conversely,  the  procedure  of  reference  2  can  be  applied  to  estimate 
these  parameters  from  information  obtained  by  loading  glass  panels  up  to  the 
failure  point. 

The  procedure  of  reference  2  requires  the  transformation  of  the  stresses 
induced  by  actual,  time-dependent  loads  into  nominal  stresses  corresponding  to 
a  1-minute  constant  load.     It  is  suggested  in  reference  2  that  the 
transformation  be  carried  out  by  using  the  following  relationship: 

J^f        (M,t)dt  l/n 
^60(M)  =  ]  (1) 

where  a5o(M)  =  equivalent  one-minute  stress  at  point  M,  a(M,t)  =  stress  induced 
by  actual  load  at  point  M  and  time  t,  tf  =  duration  of  loading  or  time  to 
failure  in  seconds,  whichever  is  smaller,  and  n  =  material  constant.  However, 
the  time  to  failure  t^  in  equation  1  is  unknown  for  panels  that  would  fail 
under  wind  loads.     To  the  extent  that  an  arbitrary  value  for  tf  is  used  in 
equation  1 ,  and  that  this  value  could  differ  by  as  much  as  one  or  even  two 
more  orders  of  magnitude  from  actual  values,  significant  errors  would  be 
introduced  in  estimating  a^Q(M).     Note  also  that  if  equation  1  were  used,  it 
would  be  necessary  to  evaluate  the  integral  at  each  point  M,  since  the  system 
is  nonlinear  (i.e.,  in  general  a(M,t)  is  not  proportional  to  the  wind  loading 
at  time  t,  p(t)). 
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On  account  of  the  difficulty — or  impossibility — of  specifying  tf,  and  of 
estimating  the  integrals  in  equation  1,  to  the  writers'  knowledge,  no  previous 
attempts  have  been  made  to  study  the  behavior  of  glass  cladding  subjected  to 
fluctuating  wind  loads.     Instead,  previous  research  has  focused  on  the  study  of 
the  behavior  of  glass  cladding  subjected  to  nominal  1-minute  loads  purported  to 
be  representative  of  fluctuating  loads. 

The  purpose  of  this  paper  is  to  present  a  procedure  for  studying  the  behavior 
of  glass  under  arbitrary  time-dependent  loads,  including  fluctuating  wind 
loads.     Like  reference  2,  the  present  procedure  is  based  on  a  nonlinear  analysis 
of  stresses  that  develop  under  the  action  of  the  external  loads,  and  on  a 
Weibull  probabilistic  model  for  the  strength  of  glass.     However,  unlike 
reference  2,  the  procedure  presented  here  incorporates  phenomenological  models 
describing  the  fracture  mechanisms  of  glass  developed  in  the  last  decade  by 
Wiederhorn  and  other  workers  (e.g.,  see  references  5,  25,  26),  as  well  as 
information  on  the  fluctuating  character  of  the  wind  loads.     No  a  priori 
assumptions  are  needed  with  regard  to  the  time  to  failure,  tf,  which  is  one  of 
the  outputs  of  the  procedure.     As  in  reference  2,  it  is  assumed  that  temperature 
and  humidity  effects  can  be  neglected.     Additional  outputs  include  the  location 
and  direction  of  the  failure  initiation  crack,  and  the  amount  of  strength 
degradation  due  to  the  action  of  the  fluctuating  wind  load. 

To  provide  a  background  for  the  development  of  the  proposed  procedure,  basic 
elements  of  the  fracture  mechanics  of  glass  will  be  briefly  summarized.  A 
method  for  obtaining  time-dependent  stresses  from  time-dependent  loads,  which 
utilizes  a  computer  program  developed  at  Texas  Tech  University,  will  be 
described.     A  brief  section  will  be  devoted  to  the  subject  of  time-dependent 
wind  loading  on  cladding  glass.     The  proposed  procedure  for  investigating  glass 
behavior  under  arbitrary  loads  will- then  be  presented.     The  procedure  will  be 
applied  to  obtain  estimates  of  probabilities  of  failure  of  a  glass  panel  under 
fluctuating  wind  loads  and  under  one-minute  constant  loads.     Such  estimates  will 
then  be  used  to  illustrate  a  method  for  assessing  current  practices  for  the 
design  of  glass  cladding  subjected  to  wind  loads. 
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2.     FRACTURE  MECHANICS  OF  GLASS 


The  basic  criterion  for  fracture  is  derived  from  the  Griffith  equilibrium 
expression,  which  can  be  written  as 

Ki  =  Kic  (2) 

where  Kj  =  stress  intensity  factor,  and  Kj^  =  critical  value  of  Kj .  If 
equation  2  holds,  the  system  reaches  the  state  of  instability  wherein  the  rate 
of  crack  growth  becomes  for  practical  purposes  infinite  [5,  7]  and  failure 
occurs.     Kjc  is  a  property  of  the  material  and  is  determined  experimentally. 
The  stress  intensity  factor,  Kj ,  is  proportional  to  the  actual  stresses  in 
the  material  in  the  presence  of  cracks  causing  stress  concentrations.     %  can 
be  expressed  as  follows  [5,  7]: 


Ki(t)  =  Ya(t)  /c(t)  (3) 

where 

Y  =  geometric  shape  factor 

a  =  nominal  stress  (i.e.,  stress  calculated  by  assuming  the  absence  of 
cracks ) 

c  =  crack  length 

t  =  time 

The  geometric  shape  factor,  Y,  is  assumed  to  be  constant;  this  is  equivalent 
to  assuming  that  the  crack  geometry  does  not  change  and  can  be  characterized 
by  one  dimension,  c.     According  to  experiments  reported  in  references  5  and 
22-26,  the  following  relationship  holds  for  the  rate  of  subcritical  crack 
growth  (figure  1): 


^  =  AK.(t) 
dt  I 


(4) 


The  parameters  A  and  n  depend  upon  ambient  humidity  and  temperature  and  are 
obtained  experimentally.     Equation  4  expresses  quantitatively  the  fact  that  the 
cracks  of  an  element  of  glass  subjected  to  stress  for  some  length  of  time 
will  grow  -  albeit  not  catastrophically  -  provided  that  the  stress  is  contained 
within  a  certain  range.     This  phenomenon  is  referred  to  as  static  or  dynamic 
fatigue  according  to  whether  the  stress  is  constant  or  time-dependent. 

It  follows  from  equations  2  and  3  that  the  strength  of  glass,  S,  i.e.,  the 
value  of  the  nominal  stress  at  which  failure  occurs,  is: 

S(t)  =-^M=^  (5) 
Y/c(t) 


3 


If  Kj  and  c  are  eliminated  from  equations  3,  4,  and  5  and  the  notation 
S(0)  =  Si  is  used  (S^  =  initial  strength),  the  following  relationship  is 


obtained: 
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S(t)  =  [Si      -i-  /  (t)  dT]  (6) 

where  1/B  =  (n-2)  AY^  [5].     It  follows  from  the  definition  of  S(t) 

that  failure  occurs  if 

a(t)  >  S(t)  (7) 

In  the  calculations  presented  in  this  paper  it  is  assumed  that  the  variability 
of  A  and  n  in  equation  4  is  small  and  that  its  effect  upon  the  results  being 
sought  can  be  neglected.     As  far  as  the  variability  of  B  is  concerned,  the 
following  comment  is  in  order.     Because  B  is  a  function  of  Y,  which  in  turn 
depends  upon  the  flaw  shape,  its  variability  is  difficult — if  at  all  possible — 
to  ascertain.     If  a  deterministic  value  of  B  is  assumed,  based  on  a  conven- 
tional value  of  Y,  it  is  possible  by  using  equation  6  to  obtain  values  of  Si 
from  experiments  in  which  n.  A,  Kj^,,  a(t)  are  known,  and  the  time  to  failure, 
tf,  and  the  strength,  S(tf),  are  measured.     The  empirical  probability  distri- 
butions of  the  initial  strength  estimated  from  the  values  S-^  so  obtained 
automatically  reflect  the  actual  shapes  of  the  flaws. 


The  initial  strength  Si  of  an  element  with  area,  a,  experiencing  a  uniform, 
direction-independent  state  of  tensile  stress  throughout  one  of  its  outer 
faces  can  be  described  probabilistically  by  a  Weibull  distribution  [21],  i.e., 

,    ,  Sj  .m, 

P(si,a)  =  1  -  exp  {-( — i—)  }  for  m  >  1  (8) 

So(a) 

where  P(si,a)  =  probability  that  the  random  variable  S^  <  s^,  So(a)  =  scale 
parameter  (characteristic  strength),  and  m  =  shape  (tail  length)  parameter. 

The  scale  parameters,  SQ(a)  and  SQ(aj),  of  two  elements  with  areas  a  and  aj, 
respectively,  each  experiencing  uniform  direction-independent  states  of  tensile 
stress  throughout  one  of  its  outer  faces,  can  be  written  as  [6]: 

J_ 
a,  m 

S^(a)  =  S^(a,)(-J^)  (9) 
o  o     i  a 


This  relation  reflects  the  dependence  of  strength  distribution  upon  the 
distribution  of  flaw  lengths.     The  larger  the  area  of  the  panel,  the  larger  will 
be  the  number  of  flaws  and,  therefore,  the  larger  the  probability  that  the  area 
will  contain  a  severe  flaw  to  which,  by  virtue  of  equation  5,  there  corresponds 
a  relatively  low  initial  strength. 

Consider  now  an  element  of  glass  with  uniform  but  direction-dependent  stresses 
throughout  one  of  its  outer  faces.     In  this  case,  shear  stresses  are  present, 


4 


in  addition  to  normal  stresses.     The  effect  of  the  normal  stresses  is  by  far 
the  strongest,  however,  as  far  as  crack  propagation  is  concerned  [7,  p.  52], 
Therefore  it  may  be  assumed  to  a  first  approximation — as  is  done  in 
reference  2 — that  the  effects  of  the  shear  stresses  can  be  neglected.     In  the 
case  of  the  element  now  being  considered  failure  will  not  necessarily  be  ini- 
tiated by  a  flaw  normal  to  the  maximum  principal  stress  -  indeed,  it  may  well 
happen  that  all  such  flaws  are  relatively  small.     Neither  will  failure  be 
necessarily  initiated  by  the  largest  flaw  within  the  element,  since  that  flaw 
may  well  be  perpendicular  to  a  relatively  low  normal  stress.     Rather,  failure 
will  be  initiated  by  the  largest  of  the  flaws  oriented  along  some  direction,  af , 
such  that 

an  (t.  Of)  >  ^  (10) 

where        (t,  af)  =  normal  stress  perpendicular  to  direction  af  (equations 

5  and  7)  and  Cj^^x         ^^f)  ~  length  of  largest  flaw  oriented  along  direction  af . 

Under  the  assumption  that  the  distribution  of  flaw  orientation  is  uniform  [2, 
p.  84],  equation  9  can  be  applied  to  the  case  of  an  element  subjected  to  uni- 
form but  direction-dependent  state  of  stress  as  follows.     The  probability  of 
occurrence  within  a  certain  area  of  flaws  with  crack  length  dimension  less 
than  c  and  having  orientation  angles  a^  <  a  <  a2  is  equal  to  (a2  ~  a^)  times 
the  probability  of  occurrence  within  an  area  of  flaws  with  dimension  less  than 
c  and  having  any  orientation  0  <  a  <2tt.     The  following  relation  is  therefore 
consistent  with  equation  9  for  any  elemental  circular  area: 


So(ai  <  a  <  a2)  =  8^(0  <  a  <  2tt)  (  H-  ) 

^2  ^1 


(11) 
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3.     STRESS  ANALYSIS 


Equations  6  and  7  clearly  show  that  to  Investigate  the  behavior  of  glass, 
information  is  required  on  the  time  history  of  the  stresses  induced  in  the 
glass  by  the  external  loads. 

The  out-of-plane  deflections  of  glass  plates  subjected  to  lateral  loads  can  be 
large  relative  to  the  thickness  of  the  plate.     The  plate  develops  substantial 
mid-plane  membrane  stresses  in  this  condition,  and  the  von  Karman 
equations  [18,  20]  must  be  used  to  account  for  this  effect.     Based  on  these 
equations,  Vallabhan  and  Wang  [19]  have  developed  a  program  that  uses  the 
finite  difference  method  to  determine  the  deflections  and  stresses  in  uniformly 
loaded  simply-supported  thin  glass  plates  having  boundary  conditions  allowing 
for  in-plane  movement.     With  the  permission  of  the  Institute  for  Disaster 
Research,  Texas  Tech  University,  this  program  was  employed  in  the  present  study 
to  determine  the  stresses  at  various  locations  on  the  plate.     The  results  from 
the  program  were  non-dimensionalized  as  follows  [9] 

LF  =  pb^/Dh  (12) 

SF  =  ab2h/D  (13) 


where  LF  and  SF  =  loading  and  stress  factors,  respectively,  p  =  uniform  pressure 
loading,   a  =  corresponding  stress,  b  =  smaller  side  of  rectangular  plate, 
D  =  flexural  rigidity,  defined  by 


12(l-v^) 

E  =  modulus  of  elasticity,  v  =  Poisson's  ratio,  and  h  =  thickness  of  plate. 

For  a  square  plate,  the  relationship  between  SF  and  LF  for  various  locations 
is  shown  in  figure  2.     Note  that  the  relationship  is  closer  to  being  linear 
for  the  corner  than  for  the  center  stresses. 

Relationships  for  each  location  on  the  plate  from  the  grid  used  for  the  finite 
difference  analysis  were  determined  and  represented  analytically  in  terms  of 
piece-wise  linear  curves.     On  the  basis  of  such  relationships,  given  a  time- 
dependent  loading,  the  plate  dimensions,  and  the  material  properties,  the 
principal  stresses  and  normal  stresses  corresponding  to  various  directions  can 
be  calculated  as  functions  of  time. 
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4.     WIND  LOADING 


Wind  loading  on  cladding  has  a  fluctuating  character.     Because  the  dimensions 
of  cladding  panels  are  relatively  small  (on  the  order  of  a  few  meters  at  most) 
it  is  acceptable  to  assume  that  the  fluctuating  loads  acting  over  the  area  of 
the  panel  are  uniformly  distributed  at  any  one  instant  and  proportional  to  the 
pressures  at  the  panel  center. 

The  wind  loads  depend  upon  the  local  extreme  wind  climate  [16,  17],  the  features 
of  the  oncoming  air  flow  (which  are  in  turn  dependent  upon  the  roughness  of 
the  surrounding  terrain  and  the  possible  presence  of  neighboring  buildings), 
and  the  aerodynamic  characteristics  of  the  building  in  question.     These  charac- 
teristics differ  from  point  to  point  on  the  building  facades  and  are,  in 
addition,  dependent  upon  wind  direction. 

Wind  loading  time  histories  are  obtained  principally  from  wind  tunnel  tests. 
Any  given  loading  time  history  can  be  analyzed  to  determine  its  statistical 
characteristics;  e.g.,  mean,  standard  deviation  and  peak.     Models  of  the  loading 
can  be  obtained  in  the  frequency  domain  (using  a  spectral  approach,  for  example, 
see  reference  17),  or  in  the  time  domain,  e.g.,  using  a  Box-Jenkins  approach 
[14],     Knowledge  of  these  models  in  turn  allows  the  numerical  simulation  of 
loading  time  histories  that  correspond  to  various  reference  wind  speeds.  In 
the  particular  case  of  a  mean  wind  speed  normal  to  a  building  face  and  in  the 
absence  of  neighboring  buildings  which  significantly  alter  the  oncoming  flow 
field,  at  a  point  Q  near  the  center  of  the  building  face  the  pressures  may  be 
written  approximately  as 


where  p  and  p^  =  mean  and  fluctuating  pressure,  respectively,  p  =  air  density, 
Cp(Q)  =  pressure  coefficient  at  point  Q,  z  =  elevation  of  point  Q,  v(z)  =  mean 
wind  speed,  at  elevation  z  and  v^(z,t)  =  longitudinal  fluctuating  component  of 
wind  speed.     The  fluctuating  pressures  can  be  simulated  from  spectral 
information  or  v'(z,t)  [17]. 


P  (Q)  =  2  P  Cp(Q)  V  (z) 


(14) 


p'(Q,t)  =  p  Cp(Q)  v(z)v'(z,t) 


(15) 
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5.     PROPOSED  PROCEDURE  FOR  INVESTIGATING  GLASS  BEHAVIOR 


Earlier  in  this  paper  it  was  noted  that  a  glass  panel  will  not  necessarily  fail 
at  the  point  of  maximum  stress  or  of  minimum  strength.     Rather,  failure  will 
occur  where  the  relationship  between  stress  and  strength  is  such  that  equation  7 
is  satisfied,  where  S(t)  is  given  by  equation  6,  and  S±  is  described  probabil- 
istically by  equation  8  (recall  that  the  effects  of  area  and  of  stress 
directionality  are  included  in  the  procedure  through  equations  9  and  11). 

With  this  background  it  is  now  possible  to  describe  the  proposed  procedure  for 
investigating  glass  behavior.     The  procedure  entails  the  following  steps: 

1.  For  any  given  mean  wind  speed  and  direction  generate  the  time  history 
of  the  wind  loading  for  the  glass  cladding  panel  of  concern  by  Monte 
Carlo  simulation  as  outlined  in  the  section  "Wind  Loading." 

2.  Using  the  procedure  outlined  in  the  section  "Stress  Analysis"  obtain 

from  the  time  history  of  the  wind  loading  the  time  histories  of  the 

stresses  normal  to  the  directions        =  — ,   (k  =  1,  2,   ...  r)  at  the 

2r 

center,       ,  of  each  of  the  elements  into  which  the  panel  is  divided 
for  numerical  computation  purposes. 

3.  From  equations  8  and  9  generate  by  Monte  Carlo  simulation  initial 
strengths  Si  (Mj »  o^k.)*  where  a]f_  =  kAa  (k  =  0,  1,  2,  . . .  n)  -  see 
figure  3. 

4.  Evaluate  numerically  the  expression 

S(M.,   o,^,  t)  =  {Si  ^(M.,  a^)  -4-  /  [o'^in.,  a^,   x)]  dj}^  ^  (16) 

B    Q  J 

for  all  j  and  k  and  t  =  mAt ,  where  At  =  incremental  time  used  in 
numerical  computations,  and  m  =  1,  2,   ..,  . 

5.  Computation  is  stopped  if  for  at  any  set  j,  k, 

a(Mj  ,   aj^,  mAt)  >_  S(Mj,   aj^ ,  mAt).  (17) 
At  that  point,  failure  of  the  general  panel  has  occurred.  ^ 
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6.     INPUT  DATA  FOR  NUMERICAL  CALCULATIONS 


Numerical  examples  were  carried  out  for  the  case  of  annealed  float  glass 
panels  simply  supported  on  four  sides  with  dimensions  4  ft  x4  ft  x  1/8  in 
(1.2  mx  1.2  mx  3  mm).     The  following  parameters  were  used:     S^(a^  =  1  m  ) 

n 

=  35.2  MPa  [4] ,  m  =  6  [4] ,  A  =  1.08  (MPa)-n  ml-2  sec"!   [23],  n  =  19.69  [23], 

Ki    =  0.75  MPa  [23],  v  =  0.22  [2] ,  E  =  0.0689  MPa  [2].     The  values  A  and  n 
C 

correspond  to  50  percent  humidity  and  20°C  temperature  (see  figure  1).     It  was 
assumed  that  the  geometric  shape  factor  Y  =  1.12  [see  reference  10  for  details. 
It  was  further  assumed  that  the  glass  panel  is  located  near  the  center  of  a 
building  facade  at  an  elevation  of  z  =  150  ft  (46  m) ,  that  it  is  subjected  to 
winds  normal  to  the  building  face,  that  the  building  is  located  in  open 
terrain,  and  that  Cp  =  1.0  in  equations  14  and  15. 

For  any  given  mean  wind  speed,  v,  the  corresponding  mean  wind  pressure  was 
obtained  by  using  equation  14.     The  fluctuating  wind  pressures  were  generated 
by  numerical  simulation  using  equation  15  and  the  following  expression  for  the 
spectrum  of  the  longitudinal  velocity  fluctuations  v',  [17]: 


n  S„  (n) 


v  _  200f 


u2  (1  +  50f)5/3 


(18) 


where  u    =  v(z)/[2.5  i!,n(z/zo)],     n  =  frequency,  f  =  nz/v(z) ,  z  =  elevation  of 
* 

panel  in  meters,  and  the  roughness  length  corresponding  to  open  terrain  Zq  = 
0.07  m  [17]. 
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7.     PROBABILITIES  OF  GLASS  FAILURE 


Failure  of  Glass  Cladding  Under  Fluctuating  Wind  Loads.     Consider  the  wind 
load,  p-i^(t),  corresponding  to  a  mean  wind  speed,  v(z).     The  subscript,  1,  indi- 
cates that  p±(t)  is  the  i-th  realization  of  the  specified  stochastic  process. 
Let  Xp,  Yp  be  the  coordinates  of  the  centers  of  the  elements  into  which  the 
glass  panel  is  divided.     The  size  of  these  elements  is  sufficiently  small  so 
that  the  variation  of  the  stresses  over  each  element  is  insignificant  for  prac- 
tical purposes.     Let  o^.  be  the  directions,  perpendicular  to  the  normal  stresses, 
being  considered  at  each  point.     Obtain  the  stresses  a-j^(xp,  yq,  a^- ,  t)  corre- 
sponding to  Pi(t)  for  all  values  p,  q,  r.     Assume  that  the  number  of  panels 
(the  sample  size)  is  L.     Simulate  initial  strengths  Sqj^^  (xp,  yq,  aj-)  for  one 
given  set  of  values  p,  q,  r  for  all  the  panels  of  the  sample  (l  =  1,  2,  ...  L). 
Rank  order  these  strengths  beginning  with  the  smallest.     Denote  this  smallest 
value  by  SqI  (xp,  yq,  a^.).     Compare        (xp,  yq,  aj-,  t)  to  S]^(t)  by  using 
equation  6.     If  a-j^  ^  Sj^  for  every  t,  then  no  panel  will  break  at  point  Xp,  yq 
in  the  direction  of  a^..     However,  if        >        for  some  t,  the  panel  for  wnich 
the  initial  strength  at  Xp,  y^,  a^.  is  S^i  will  break.     One  breakage  is  recorded, 
and  the  panel  concerned  is  eliminated  from  further  consideration.     If  a±  >  S2, 
where  S2  is  derived  from  the  second  smallest  initial  strength  for  the  given 
p,  q,  r,  then  breakage  initiated  at  Xp,  yq,   a^-  will  occur.     Repeat  the  compari- 
son until  <  Sj^.     The  number  of  failures  initiated  at  Xp,  yq,  a^-  will  then 

be  N=N-1 ,  and  the  corresponding  panels  are  eliminated  from  further  considera- 
tion.    Repeat  the  procedure  for  all  values,  p,  q,  r.     The  probability  of  failure 

under  load  Pj:(t)  is  F[f allure  |  pj^(t)]  =    E(N/L)j^.     Repeat  the  the  procedure 
for  different  realizations  of  the  pressure.     The  probability  of  failure  under 

G 

load  p(t)  is  P[fallure   |  p(t)]  =  E  Prob.   [failure   |  pi(t)]/G,  where  the 

1 

summation  is  over  the  number  of  realizations,  G,  being  considered. 

A  cumulative  distribution  function  P[failure   |   p(t)]  is  derived  indicating 
that  the  strength  of  the  plate  is  less  than  required  to  withstand  the  load  p(t) 
in  100  X  P[fallure  |  p(t)]  percent  of  the  cases.     Note  that  P[f allure  |  p(t)] 
embodies  both  the  characteristics  of  the  material  and  the  aerodynamic  character- 
istics of  the  loading  under  consideration.     These  two  types  of  characteristics 
are  inseparable  owing  to  the  dependence  of  glass  strength  upon  load  time 
history. 

For  convenience,  the  load  p(t)  may  be  indexed  by  its  mean  value  p,  so  that  the 
notation  P[failure  |   p],  or  simply  P^  (p)  may  be  substituted  for  P[fallure 
I  p(t)].     Estimates  of  points  of  the  cumulative  distribution  function  Pf(p) 
estimated  for  the  conditions  described  in  the  section  "Input  Data  for  Numerical 
Calculations"  are  shown  in  figure  4.     These  points  were  obtained  for  simulated 
failure  tests  carried  out  on  two  sets  of  1,000  panels,  each  set  being  subjected 
to  different  realizations  of  the  fluctuating  pressure,  p(t). 

Typical  time  histories  of  pressure,  stress,  and  strength  at  panel  points  where 
failure  was  found  to  occur  are  shown  in  figure  5.     Note  the  continuous  strength 
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degradation  under  load.     In  certain  instances  a  storm  can  cause  significant 
strength  degradation  without  causing  failure.     The  weakened  panel  could  then 
break  under  the  action  of  subsequent,  less  intense  storms,  should  such  storms 
occur  during  its  lifetime. 

Failure  of  Glass  Cladding  Under  1-Minute  Constant  Load.     Let  the  normal  stress 
induced  by  the  1-minute  load,  P60>  the  initial  strength  be  denoted  by 

a5o(M,  a)  and  S-j^(M,  a),  respectively,  where  M  and  a  are  the  point  and  the 
direction  under  consideration.     It  is  assumed  that  the  60-second  action  of  059 
(M,   a)  causes  failure.     It  follov7s  from  equation  6  that,  with  negligible  error, 
the  following  relation  holds: 

S.     (M,   a)  ^/'^ 


The  procedure  for  obtaining  pgQ  is  as  follows:     Generate  initial  strengths, 
^ol  (^p»  yq>  ^'r)  ^  ~  ^'  2,...  L  panels.     Use  equation  19  to  solve  for  059 

(M,  a)  at  each  Xp,  yq,  a^.  for  panel  i.     Calculate  the  value  of  P50  correspon- 
ding to  each  a^QvM,   a).     Find  the  minimum  pgQ  foi"  panel  i;  it  represents  the 
lowest  60-second  pressure  loading  for  which  failure  of  panel  i  will  occur. 
Repeat  the  procedure  for  i  =  I,  2,   ...  L  panels.     The  results  obtained  can  be 
plotted  in  the  form  of  a  cumulative  distribution  function  P [failure  |  P6oK 
denoted  as  Pf(p5o)«     Figure  6  shows  a  cumulative  distribution  P^Cp^q)  for  the 
conditions  described  in  the  section  "Input  Data  for  Numerical  Calculations," 
based  results  of  obtained  for  a  set  of  4,000  panels. 

Note  that  according  to  the  PPG  charts   [2]  the  1-minute  load,  P60»  corresponding 
to  a  probability  of  failure  of  8  in  1,000  is  23  psf.     The  corresponding  value 
of  p5o  indicated  by  figure  6  is  about  19  psf. 
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8.     ASSESSMENT  OF  CURRENT  PROCEDURES  FOR  DESIGNING  GLASS  CLADDING 


The  purpose  of  this  section  is  to  compare  the  current  design  loadings  with  the 
present  procedure. 

Procedure  Used  in  Refs.  12-13.     This  procedure  is  based  on  the  assumption  that 
the  effect  of  wind  loading  is  determined  solely  by  the  peak  load  averaged  over 
1  or  2  seconds  that  occurs  during  the  storm.     It  is  assumed  that  the  loading 
on  the  cladding  which  is  subjected  to  a  fluctuating  wind  pressure  is  equivalent 

C 

to  a  1-minute  constant  pressure,  P50>  defined  as 

J_ 

P60  =  Ppk  (^>'  (20) 

where  pp]^  =  peak  pressure  averaged  over  the  time,  tp^,  that  occurs  during  the 
storm  being  considered,  and  tp]^  =  1-2  seconds.     According  to  this  assumption, 
then,  provided  that  Pp^  is  the  same,  it  does  not  matter  whether  the  mean 
loading  is  large  and  the  fluctuations  are  small  or  vice-versa.     The  probability 
of  failure  of  a  panel  subjected  to  a  load,  p(t),  implicit  in  reference  12  is 
thus  equal  to  the  probability  of  failure  of  that  panel  under  the  action  of  a 
1-minute  load,  P^q,  obtained  from  p(t)  by  using  the  above  equation. 

Procedure  Proposed  in  Ref.  3.     This  procedure  assumes  that  the  equivalent 
1-minute  pressure  loading  is  given  by 

t 

D  r  ^  n  1/n 

P60  =  r  P  eo''  ^''^ 

where  tg  =  duration  of  storm. 

Comparison  of  Design  Loadings  Based  on  Various  Procedures.     The  purpose  of 

C  D 

this  section  is  to  compare  the  loads  p(t)  ,  P50»  and  p^Q  corresponding  to  various 
probabilities  of  failure.     These  loads  are  obtained  from  the  cumulative 
distribution  functions  of  figs.  4  and  6. 

To  the  storm  with  mean  speed  v  these  corresponds  a  mean  pressure  p,  a 

C 

1-minute  load  obtained  in  accordance  with  reference  12,  PgQ*  ^  1-minute  load 

D 

obtained  in  accordance  with  reference  3,  PgQ*  conditions  described  in 

C    _  D  _ 

the  section  "Input  Data  for  Numerical  Calculations,"  the  ratios  P50/P  P60^P 

C 

were  found    to  be  approximately  1.08  and  1.26,  respectively.     Values  of  pgQ 

D  _ 
and  pgQ  corresponding  to  various  values  p,  and  the  probabilities  of  failure 
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C  D 
Pf(p)   =  P[failure   |  p(t)],  Pf  (p^g  =  P60)'  ^^^^        (P60      P60)'  obtained  from 
figs.  4  and  6,  are  listed  in  table  1.     It  is  seen  that,  in  this  instance,  the 
estimates  based  on  the  procedure  of  reference  12  appear  to  be  overly  optimistic, 
i.e.,  they  appear  to  underestimate  the  probability  of  failure  of  the  panel 
under  any  given  storm  and,  therefore,  the  probability  of  failure  of  the  panel 
during  its  lifetime.     The  probability  estimates  based  on  the  procedure  of 
reference  3  are  somewhat  closer  to  those  based  on  the  time  history  of  the 
stresses.     Note  that  for  storms  causing  rates  of  failure  of  about  8  in  a  1,000 
the  probability  estimates  based  on  the  stress  time  history  and  on  the  nominal 
C  D 

1-minute  loads  p^g  ^^'^  P60  happen  in  this  case  to  be  relatively  close.  The 
respective  discrepancies  increase  considerably  in  the  case  of  stronger  storms. 
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9.  CONCLUSIONS 


A  procedure  for  investigating  glass  cladding  behavior  under  arbitrary  loads, 
including  fluctuating  wind  loads,  was  presented.     The  procedure  accounts  for 
the  fact  that  internal  stresses  are  nonlinear  functions  of  the  external  loads, 
that  initial  glass  strengths  are  random  functions  of  position  and  direction, 
and  that  the  glass  strength  undergoes  degradation  under  the  action  of  external 
loads  in  accordance  with  basic  fracture  mechanics  laws  that  reflect  subcritical 
crack  growth.     Numerical  examples  were  presented,  and  corresponding  probability 
distributions  were  calculated,  indicating  the  probability  of  failure  of  a 
specified  panel  subjected  to  fluctuating  wind  loads  and  to  1-minute  constant 
loads.     These  curves  are  used  to  illustrate  a  method  for  assessing  current 
glass  cladding  design  procedures.     For  the  case  considered  in  the  paper  it  was 
found  that  procedures  based  on  the  transformation  of  the  wind  load  averaged 
over  1-2  seconds  into  an  equivalent  1-minute  load  appear  to  result  in  overly 
optimistic  assessments  of  the  probability  of  failure  of  glass  cladding  under 
wind  loads.     The  work  reported  in  the  paper  is  part  of  an  ongoing  window 
cladding  research  program  being  conducted  at  the  National  Bureau  of  Standards. 
Future  work  will  consider  the  effect  of  fluctuating  loads  on  corner  and  eave 
panels,  and  the  effect  upon  fracture  load  predictions  of  the  variability  of 
the  parameters  that  control  glass  behavior  under  load. 
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APPENDIX  I.  NOTATION 


a,  aj  =  area  under  uniform  tension 

A  =  material  constant  relating  rate  of  crack  propagation  to  the  stress 
intensity  factor  Kj 

B  =  material  constant  used  in  fracture  mechanics  analysis 

b  =  length  of  side 

c  =  flaw  (crack)  length 

D  =  flexural  rigidity 

E  =  modulus  of  elasticity 

f  =  Monin  or  similarity  coordinate 

G  =  total  number  of  storm  realizations 

h  =  plate  thickness 

Kj  =  stress  intensity  factor 

Kj  =  critical  value  of  Kj 
c 

LF  =  load  factor 

m  =  Weibull  parameter 

n  =  material  constant 

n  =  frequency  in  Hz 

N  =  total  number  of  panel  failures /storm  realization 

p(t)  =  wind  pressure  loading  at  time  t 

p(Q)  =  mean  component  of  the  quasi-static  pressure  at  point  Q 

p(Q,t)  =  fluctuating  component  of  the  quasi-static  pressure  at  point  Q 

P50  ~  orie  minute  effective  loading 
C 

P50  ~  oi^G  minute  effective  loading  defined  by  equation  20 
D 

P50  =  one  minute  effective  loading  defined  by  equation  21 
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Pplj^  =  peak  pressure 

=  Internal  velocity  pressure 

qp  =  effective  velocity  pressure  for  parts  and  portions 

s  =  sample  standard  deviation 

S(t)  =  strength  at  time  t 

Si  =  initial  strength 

So  =  Weibull  location  parameter 

SF  =  stress  factor 

t  =  time 

tpk  =  time  duration  of  peak  pressure 

Vf  =  fastest  mile  wind  speed 

X  =  sample  mean 

Y  =  proportionality  constant 

y  =  parameter  calculated  in  Type  I  distribution 

z  =  elevation 

Of  =  direction 

V  =  Poisson's  ratio 
p  =  density  of  air 

a(M,t)  =  nominal  stress  induced  by  a  pressure  load  at  point  M  and  time  t 
Oj^Ct,   Of)  =  normal  stress  perpendicular  to  direction  af  at  time  t 

%0^^)  ~  equivalent  one-minute  stress  at  point  M 
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APPENDIX  II.     COMPUTER  LISTING 


c 

C  GLASS  ANALYSIS  PROGRAM 
C 

c 

C  THIS  PROGRAM  ESTIMATES  THE  PROBABILITY  OF  FAILURE 

C  OF  SQUARE  GLASS  PANELS  SUBJECTED  TO  ARBITRARY 

C  PRESSURE  LOADINGS.     THE  PROGRAM  CAN  BE  EASILY 

C  MODIFIED  TO  ACCOMODATE  RECTANGULAR  PANELS. 

C  THE  LOADINGS  CAN  BE  CONSTANT  ONE-MINUTE  LOADS, 

C  DETERMINISTIC  TIME  DEPENDENT  LOADS,  OR 

C  FLUCTUATING  LOADS.     THE  LATTER  ARE  SIMULATED 

C  USING  QUASI-STATIC  THEORY  AND  AUTOREGRESSIVE  MODELS, 

C  OR  CAN  BE  READ  IN  AS  FLUCTUATING  OR  CONSTANT. 

C  THE  INITIAL  STRENGTH  VALUES  FOR  EACH  LOCATION 

C  ON  THE  PANEL  ARE  GENERATED  FROM  INPUT  WEIBULL 

C  PARAMETERS,  OR  CAN  BE  READ  IN,  OR  SIMPLY  INPUT 

C  AS  CONSTANT.     STRENGTH  DEGRADATION  IS  DETERMINED 

C  AND  IF  DESIRED,  THE  STRESS  AND  STRENGTH  TIME 

C  HISTORIES  OF  THOSE  PANELS  WITH  DEGRADATION 

C  OF  A  CERTAIN  PERCENTAGE  CAN  BE  WRITTEN  ONTO  A 

C  FILE. 

C 

C 

C  THE  SIXTY-SECOND  LOADINGS  GIVEN  BY  EQNS.  20  AND  21  ARE 

C  CALCULATED  FOR  EACH  P(T)   (FLUCTUATING  LOADING). 

C 

C 

c 

C  TO  OBTAIN  FAILURE  STATISTICS,  NUMERICAL  EXPERIMENTS 

G  ARE  CONDUCTED  ... 

C  UP  TO  1000  PANELS  CAN  BE  SUBJECTED  TO  AS  MANY 

C  AS  10  STORMS,  EACH  HAVING  THE  SAME  MEAN  BUT 

C  WITH  FLUCTUATINF  TIME  HISTORIES  THAT  ARE 

C  DIFFERENT  ON  ACCOUNT  OF  THE  RANDOM  NATURE  OF 

C  THE  FLUCTUATIONS.     EACH  STORM  CAN  BE  OF  DURATION 

C  900  SECONDS  OR  LESS.     TO  SUMMARIZE... 

C  TEN  STORMS  CAN  BE  TESTED. 

C  1000  PANELS/ STORM  CAN  BE  TESTED. 

C  STORM  DURATION  CAN  BE  900  SECONDS  OR  LESS. 

C  THESE  PARAMETERS  CAN  BE  CHANGED  BY  INCREASING 

C  OR  DECREASING  THE  DIMENSION  STATEMENTS  BELOW... 

C  INPUT  PARAMETERS  ARE  DEFINED  BELOW  

C 

c 

C  THE  STATEMENTS  CONCERNING  THE  PRESSURE  TIME 

C  HISTORY  IN  THIS  PROGRAM  CORRESPOND  TO  THE 

C  CASE  WHERE  ( 1 )  THE  VELOCITY  IS  PERPENDICULAR 

C  TO  THE  PLANE  OF  THE  FACADE  AND  (2)  THE  PANEL 

C  IS  SUFFICIENTLY  FAR  FROM  THE  EAVE  OF  THE 

C  GROUND  (SEE  EQNS  1M  AND  15  IN  THE  REPORT). 

C  FOR  ANY  OTHER  POSITION  OF  THE  PANEL,  THE  PRESSURES 

C  ARE  DETERMINED  IN  ACCORDANCE  WITH  THE  RELATION: 

C 

C  P(Y,Z,T)  =  0.5  RHO  (U(H)»«2  )  CP(Y,Z,T) 
C 

C  WHERE  P(Y,Z,T)  =  TIME  DEPENDENT  PRESSURE  AT  ELEVATION  Z 

C  AT  DISTANCE  Y  FROM  THE  CENTERLINE  OF 

C  THE  FACADE  AT  TIME  T 

C  RHO           =  AIR  DENSITY 

C  U(H)          =  MEAN  WIND  SPEED  AT  THE  TOP  OF  THE 

C  BUILDING 
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C  CP(Y,Z,T)=  PRESSURE  COEFFICIENT  CORRESPONDING  TO 

C  THE  WIND  TUNNEL  OR  FULL  SCALE  REFERENCED 

C  WITH  RESPECT  TO  U(H) 
C 

c 
c 
c 

c 

C  PARAMETER  DEFINITIONS  IN  ORDER  OF  INPUT... 


C 

c 

C  NSTORM  NUMBER  OF  STORMS 

C  NPANEL  NUMBER  OF  PANELS 

C  NDUR  STORM  DURATION  (  SECONDS  ) 

C  DT  TIME  STEP  (  SECONDS  ) 

C  SIDEA  LENGTH  OF  PANEL  SIDE  (INCHES) 

C  TH  THICKNESS  (INCHES) 

C  E  MODULUS  OF  ELASTICITY  (PSI) 

C  PR  POISSONS  RATIO 

C  KIC  STRESS  INTENSITY  FACTOR  (COMPATIBLE  W/  MPA  UNITS) 

C  A  FRACTURE  MECHANICS  PARAMETER (SAME  AS  WITH  KIC) 

C  RY  GEOMETRIC  SHAPE  PARAMETER 

C  RN  POWER  USED  IN  EQN.  6 

C  SO  UNITS  =  MPA,  WEIBULL  DISTRIBUTION  PARAMETER 

C  CO  Mrl/CG  IN  WEIBULL  DISTRIBUTION,  EQN,  8  IN 
C  THE  REPORT 

C  SCONT  CONSTANT  STRENGTH  VALUE  ,UNITS=MPA 

C  STCONT  CONSTANT  STRESS  VALUE,  UNITS  =  MPA 

C  DAF(1)  NAME  OF  FILE  WHICH  CONTAINS  INITIAL 

C  STRENGTH  VALUES,  IF  THIS  OPTION  IS 

C  CHOSEN 

C  ILOAD  LOAD  TYPE 

C  =  1   ,  SIMULATE  FROM  AR(2)  MODEL 

C  =  2  ,  CONSTANT,  READ  IN  FROM  LOGICAL  UNIT  5 

C  =  3  ,  READ  IN  FROM  LOGICAL  UNIT  10 

C  ISTREN  INITIAL  STRENGTH  TYPE 

C  =  1   ,  WEIBULL  DISTRIBUTION 

C  =  2  ,  CONSTANT,  INPUT  VALUE 

C  =  3  ,  READ  IN  FROM  LOGICAL  UNIT  11 

C  ISTRES  STRESS  TIME  HISTORY  TYPE 

C  =  1   ,  CALCULATE  FROM  P(T) 

C  =  2  ,  CONSTANT,  INPUT  VALUE 

C  IWRIT  PRINTING  OPTION 

C  =  1   .PRINT  ONLY  FINAL  TABLE 

C  =2  .PRINT  OUT  FAILURE  TIME  &  INITIAL  STRENGTH 
C     ■  AT  LOCATIONS  WHERE  FAILURE  OCCURS 

0  &  THE  FINAL  TABLE 

C  =3  .PRINT  OUT  FAILURE  TIME  &  INITIAL  STRENGTH 
C  AT  LOCATIONS  WHERE  FAILURE  OCCURS, 

C  AS  WELL  AS  THE  STRENGTH  &  STRESS  TIME 

C  HISTORIES,  &  THE  FINAL  TABLE 

C  IRAT  RATIO  CHECK 

C  =  1   .  CHECK  FOR  %  DEGRADATION 

C  &  PRINT  OUT  %  AT  LOCATION  &  DIRECTION 

C  =  2  .DO  NOT  CHECK 

C  =3  .CHECK  AND  PRINT  OUT  THE  STRENGTH 

C  &  STRESS  TIME  HISTORIES 

C  =  1  OR  3,  INPUT  RAT1.RAT2  %S  BELOW 

C  NOTE  RATI  MUST  BE  .GE.  RAT2 

C 

C 

c 
c 


DOUBLE  PRECISION  DSEED,TERM1 .TERM2,RDIFF.RDIR .RATIO. SLOPE 
DOUBLE  PRECISION  A1 ( 1 000) . A2( 6 ) , S( 900) , AI (6 , 1 000) ,RSIG( 36 . 6 ) 
REAL  IT,PV(10),F(1000),PMEAN,PM(10),NFAIL2(216).FMAX. 
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•UMEAN ( 1 0 ) , PF ( 1 0 ) , PCER (10), PDAG ( 1 0 ) , RAT 1 , RAT2 
DOUBLE  PRECISION  P( 900 ), SIGMA (36, 900) 
EQUIVALENCE(  P(1),S(1)  ) 

DOUBLE  PRECISION  M, COEF,D, SO , BP, A, KIC, TH, PR ,E, ROOT, RNM, RN 

DOUBLE  PRECISION  P60( 10) ,SCALE,FLEX,SC0NT,STC0N,SIGG2, 
»SINT(1000),SU,Ry 

INTEGER  MRI , ISTREN , NCOUNT , NPANEL , NDUR , I 1 , 
•ISEED , TIME ( 1 000 ) , IWRIT , IRAT , K? , K6  ( 36 ) 

CHARACTER*16  DAF(1») 

C 

C  ENTER  PARAMETERS 

C 

C 

C 

READ ( 5 , • ) NSTORM , N  PANEL , NDUR , DT 

READ(5,*)SIDEA,TH 

READ(5,»)E 

READ(5,»)PR,KIC 

READ(5,*)A,RN,RY 

READ ( 5 ,*) ILOAD , ISTREN , ISTRES , IWRIT , IRAT 
IF(  ISTREN  .EQ.   1 )READ(5,»)SU,S0,C0 
IF(  ISTREN  .EQ.  2)READ(5,*)SC0NT 
IF(  ISTRES  .EQ.  2)READ(5,*)STC0N 
IF(  ISTREN  .EQ.  3  )READ(5,8)DAF( 1 ) 
IF(  IRAT  .NE.  2  )READ(5,»)RAT1 ,RAT2 

C 

CALL  TWRITE ( ILOAD , ISTRES , ISTREN , NDUR ) 

C 

C        DEFINE  POWERS  FOR  EQUATIONS 
C 

RNN=RN+1 . 
RNM=RN-2. 
R00T=1./RNM 
R00T2=1./RN 

C 

C  CALCULATIONS 
C 

BP=(RNM»A»RY»RY*(KIC*»RNM))/2. 

FLEX= ( E  »TH»TH  »TH ) / ( 1 2 . • ( 1 - PR  »PR ) ) 

COEF=DT«BP/RNN 

SCALE=  ( SIDEA«»i| . )  /  ( FLEX'TH ) 

C 

c 
c 

CALL  GWRITE ( SIDEA , TH , E , PR , A , KIC , NPANEL , RN , BP) 

,0 

DO  5555KK=1 , NSTORM 
IF(  IWRIT  .NE.  1  )WRITE(6,7)KK 
IF(  IRAT  .NE.  2  )WRITE(9,7)KK 
IF  (  ILOAD  .EQ.  2  )  GO  TO  1111 
IF  (  ILOAD  .EQ.  3  )  GO  TO  1112 
IF(  ISTRES  .EQ.  2  )  GO  TO  2222 
IF(  ISTREN  .GT.  H  )WRITE(6 ,9000) 
IF(  ISTRES  .GT.  k  )WRITE(6 , 9000) 


c 
c 

c 
c 
c 

PARAMETER 

DEFINITIONS  FOR  P(T)  SIMULATION 

c 

CP 

STATIC  PRESSURE  COEFFICIENT 

c 

Z 

ELEVATION 

c 

ZO 

FRICTION  LENGTH 

c 

UMEAN 

MEAN  WIND  VELOCITY 

c 

MRI 

MEAN  RECURRENCE  INTERVAL 

c 
c 
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READ(5,»)CP,Z,Z0,IT,MRI 

DATA  XMIN,XMAX,NDIV/0, 100, 1000/ 

UMEAN(KK)=IT 

IT=IT«1.i»7 

C 

C    ICLOCK  IS  A  FUNCTION  WHICH  RETURNS  THE  NUMBER  OF  SECONDS 

C    SINCE  MIDNIGHT   

C 

c 

CALL  ICLOCK ( 2, ISEED) 
DSEED=DFLGAT ( ISEED ) 

C 
C 

C         GENERATE  LOAD  AND  NON-DIMENSIONALIZE 

C 

C 

CALL  PRES ( XMIN , XMAX , NDI V , DSEED ,DT,NDUR,Z,ZO,IT,CP, PMEAN , PVA , P ) 
PV(KK)  =  1i»4.«PVA 
PM(KK)  =  141|.*PMEAN 

CALL  PWRITE ( IT , MRI , DT , CP , PMEAN , P , NDUR ) 

C 
C 

c 

GO  TO 

1111  CONTINUE 
READ(5,»)P60(KK) 
DO  3331=1, NDUR 

333  P(I)=P60(KK) 
GO  TO  HkH 

1112  READ(10,»)(  P(I),  1=1   ,NDUR  ) 
C 

C 

C  FIND  THE  EQUIVALENT  SIXTY-SECOND  LOADINGS  FOR  P(T) 

C 

C 

kkh  CONTINUE 

CALL  PEQ(NDUR,DT,RN,P,P1,P2) 
PCER(KK)  =  1i<4.»P1 
PDAG(KK)=141».»P2 

C 
C 

C  CALCULATE  THE  LOAD  FACTOR  FROM  P(T)... 

C 

C 

DO  301 11=1, NDUR 
3011  P(I)=P(I)»SCALE 

C 

C  FIND  DIRECTIONAL  MULTIPLIERS  AND  CALCULATE  THE  STRESS  FACTORS... 
C 

CALL  DIR(RSIG) 

CALL  STRESS( FLEX , SIDEA , TH , NDUR , P , SIGMA ) 

C 
C 

c 

C  INITIAL  STRENGTH  DISTRIBUTION 

C  WEIBULL  DISTRIBUTION  ASSUMED  

C 

2222  CONTINUE 

IF  (  ISTRES  .NE.  2  .AND.  ISTREN  .NE.  2  )G0  TO  2223 

J=1 

JK=0 

111=0 

GO  TO  60 

2223  IF(  ILOAD  .NE.  2  .AND.  IWRIT  .NE.   1  )WRITE(6,3) 

IF(  ISTREN  .EQ.  3  )OPEN( 1 1 , IOSTAT=ISTAT, FILE=DAF( 1 ) ) 

C 
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c 

DO  111J=1,36 
111  K6(J)=0 

DO  10J=1,36 

IFdSTREN  .EQ.  2)  GO  TO  50 
IFdSTREN  .EQ.  3)  GO  TO  HO 
M=1./C0 

C 
C 

CALL  ICL0CK(2,ISEED) 

DSEED=DFLOAT(ISEED) 

CALL  GGUBS (DSEED.N PANEL, F) 

C 
C 

DO  22K=1,6 

DO  3011=1 ,NPANEL 

30  A1(I1)=SU  +  S0«((-AL0G(F(I1)))»»M) 
C 

C  RANK  THE  INITIAL  STRENGTH  VALUES  ... 
C 

C  ELIMINATE  THOSE  WHICH  ARE  TOO  LARGE  FOR  FAILURE  TO 

C  OCCUR... 

C 

C 

CALL  VSRTAD(A1,NPANEL) 

C 
C 

DO  3111=1, NPANEL 

31  AI(K,I1)=A1{I1) 

22  CONTINUE 
DO  23K=1,6 
A2(K)=AI(K,1) 

23  CONTINUE 
C 

C  VSRTAD  IS  AN  IMSL  SUPPLIED  ROUTINE  WHICH 

C    RANKS  THE  VALUES  OF  THE  INPUT  VECTOR  FROM  SMALLEST 

C    TO  LARGEST  

C 
C 

C  RANK  THE  INITIAL  STRENGTH  VALUES  ACCORDING  TO  DIRECTION. 
C 

CALL  VSRTAD (A2, 6) 

C 
C 

C  CHECK  TO  SEE  HOW  MANY  DIRECTIONS  CAN  BE  ELIMINATED  FROM 

C  COMPARISON  BELOW... 

C 

C 

DO  24K=1,6 

IF(  A2(1)   .EQ.  AI(K,1)  )K5=K 
K6(J)=K5 

24  CONTINUE 
GO  TO  10 

C 

C  READ  IN  INITIAL  STRENGTH  VALUES  AND  RANK  AS  ABOVE  FOR 

C  GENERATED  VALUES   

C 
C 
C 

40  CONTINUE 

DO  26K=1,6 

READ(11,44)   (  A1(I1),  11=1, NPANEL  ) 

C 

C  RANK  THE  INPUT  INITIAL  STRENGTH  VALUES  ... 
C 
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c 

CALL  VSRTAD(A1,NPANEL) 

C 
C 
C 

DO  2711=1 ,NPANEL 
AI(K,I1)=A1(I1) 

27  CONTINUE 
26  CONTINUE 

DO  28K=1,6 
A2(K)=AI(K,1) 

28  CONTINUE 
C 

C 

C  RANK  ACCORDING  TO  DIRECTION  

C 
C 

CALL  VSRTAD(A2,6) 

C 

C  CHECK  TO  SEE  HOW  MANY  DIRECTIONS  MAY  BE  CRITICAL... 

C 

C 

C 

DO  29K=1,6 

IF(A2(1)  .EQ.  AKK.D)  K5=K 
K6(J)=K5 

29  CONTINUE 
10  CONTINUE 
50  CONTINUE 
C 

C 

C  THE  ICQ  LOOP  IS  FOR  EACH  NPANEL  CHECK 

C 

C 

JK=0 
NPOS=36 

DO  110J=1,NPOS 

K7=K6(J) 

NCOUNT=0 

DO  20K=1,K7 

NCOUNT2=0 

JK=JK+1 

DO  100  11=1, NPANEL 

IFdSTREN  .EQ.  2)  GO  TO  60 

S(1)=A1(I1) 

GO  TO  70 
60  S(1)=SC0NT 

NCOUNT=0 

JK=JK+1 

111=111+1 
70  CONTINUE 
C 

C    EQUATION  6  IN  REPORT 

C 

C 

DO  200I=2,NDUR 
IF(  ISTRES  .EQ.  2  )  D=1.0 
IF(  ISTRES  .EQ.  2  )G0  TO  700 
SLOPE=SIGMA ( J , I ) -SIGMA ( J , I- 1 ) 
IF(  SLOPE  .EQ.  0  )  GO  TO  700 
IF(  RSIG(J,K)   .LT.  0  )G0  TO  20 
RDIR=  (RSIG(J,K))»»RN 

D=RDIR»((SIGMA(J,I)»«RNN)  -(SIGMA(J,I-1 )«»RNN)  )/SL0PE 
SIGG2=  SIGMA ( J , I ) "RSIG ( J , K ) 
GO  TO  800 

700      IF(  ISTRES  .EQ.  2  )SIGG2=STC0N 
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IF(  ISTRES  .NE.  2  )SIGG2=SIGMA(J,I)»RSIG(J,K) 

IF(  ILOAD  .EQ.   2  )D=1.0 

IF(  SIGG2  .LE.  0  )  SIGG2=0 

IF(  SIGG2  .EQ.  0  )D=0 

IF(  D  .EQ.  0  )  GO  TO  800 

D=RNN»(  SIGG2»»RN  ) 
800      TERM1=S(I-1 )»»RNM 

TERM2=C0EF»D 
C         WRITE(6,»)'  TERM1  ',TERM1 
C         WRITE(6,«)«  TERM2  ',TERM2 

RDIFF=TERM1  -TERM2 

IF(  RDIFF  .LT.  0  )  GO  TO  900 

S(I)=RDIFF»»ROOT 
200      IF(     SIGG2     .GE.  S(I)   )  GO  TO  900 
C 
C 

C  FAILURE  DOES  NOT  OCCUR  IN  THE  GIVEN  TIME  PERIOD 

C  CHECK  FOR  DEGRADATION   

C 

IF(  IRAT  .EQ.  2  )G0  TO  102 
RATIO  =  (S(NDUR)/S(1))»100. 
IF(  RATIO  .LT.  RATI   )  GO  TO  901 
102  CONTINUE 

IF(  ISTRES  .EQ.  2  .AND.  ISTREN  .EQ.  2  )  GO  TO  5555 
IF(  ISTREN  .EQ.  2  )  GO  TO  10 
GO  TO  902 

900  NCOUNT  =  NCOUNT  +  1 
NC0UNT2=NC0UNT2+1 
TIME (NCOUNT )=  I 
SINT(NCOUNT)=  S(1) 

IF(  IWRIT  .EQ.   1   .OR.  IWRIT  .EQ.  2  )  GO  TO  333^ 
NEND=I 

IF(  ISTRES  .NE.   2)  WRITE(6 , 6000) (  SIGMA(J,I) ,1=1 ,NEND) 

WRITE(6,6000)(  S(I),  I=1,NEND  ) 

IF(I1   .LT.  NPANEDGO  TO  100 

IFdSTREN  .NE.  2)  GO  TO  3334 

DO  l»Hi(3I=1  ,NEND 
HHtiS    SIGMA(1,I)=  SIGG2 

WRITE(6,6000)(  SIGMA(J,I),  I=1,NEND  ) 
3334    DO  99991= 1,NEND 
9999  S(I)=0.0 

GO  TO  100 

901  WRITE (9,*)'     RATIO  OF  FINAL  TO  INITIAL  S  =     '.RATIO,'  %  ' 
WRITE(9,»)'  AT  POSITION  '  ,  J, '  AND  DIRECTION  ',K 

IF(  ISTRES  .EQ.  2  .AND.  ISTREN  .EQ.  2  )WRITE(9,»)'  PANEL  #  ',111 

IF(ISTRES.NE.2  .AND.  ISTREN. NE. 2)WRITE(9, •) '  PANEL  #  ',11 

WRITE (9,*)'  %  CHECK  VALUES  ARE  ',RAT1,'  AND  ',RAT2 

IF  (  RATIO  .GE.  RAT2  )G0  TO  101 

IF(  IRAT  .EQ.   1   .OR,  IRAT  .EQ.  2  )G0  TO  101 

IF(  RATIO  .LT.  RAT2  )WRITE{9 ,6000) (  S(I ) , 1= 1 ,NDUR  ) 

DO  88761=1, NDUR 

8876  SIGMA(J,I)=SIGMA(J,I)»RSIG(J,K) 

IF(  RATIO  .LT.  RATI   )WRITE(9 ,6000) (  SIGMA(J,I) ,1=1 ,NDUR) 
DO  67881=1, NDUR 

6788  SIGMA(J,I)=SIGMA(J,I)/RSIG(J,K) 

101  CONTINUE 

DO  99981=1, NDUR 

9998  S(I)=0. 

C 

C 

C  FINISH  OF  100  LOOP 

C 

C 

100  CONTINUE 

C 

C 
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902      DO  99881= I.NDUR 
9988  S(I)=0. 

NFAIL2 ( JK ) =REAL ( NC0UNT2 ) 
IF(  NCOUNT  .EQ.  0  )  GO  TO  20 
IF(  IWRIT  .EQ.   1)  GO  TO  20 
WRITE(6,3) 

DO  700 11=1, NCOUNT 

7001  WRITE(6,2)I,J,K,TIME(I),SINT(I) 

IF(  ISTRES  .EQ.  2  .AND.  ISTREN  .EQ.  2)  GO  TO  5555 

C 

C  FINISH  THE  DIRECTION  LOOP 
C 

20  CONTINUE 
C 

C  FINISH  THE  LOCATION  /  POSITION  LOOP 

110  CONTINUE 

C 

C 

C  CALCULATE  THE  NUMBER  OF  FAILURES  PER  STORM  AND 

C  WRITE  TO  AN  OUTPUT  DEVICE   

C 

FMAX=0.0 
II=JK 

DO  70021=1,11 

7002  FMAX=AMAX1(FMAX,NFAIL2(I)) 
PF(KK)=FMAX/NPANEL 

IF(  ILOAD  .EQ.  2  )  GO  TO  5555 

C 

C 

C  GO  TO  NEXT  STORM... 

C 

C 

CLOSE(II) 
5555  CONTINUE 
C 

IF(  ISTREN  .EQ.  2  .AND.  ISTRES  .EQ.  2  )  GO  TO  5561 
IF(  ILOAD  .EQ.  2  )G0  TO  5558 
WRITE(6,6) 

C 

C  PRINT  THE  FINAL  PROBABILITY  ESTIMATES 

C         DEFINITIONS. . . 

C         UMEAN(K)      MEAN  WIND  SPEEDS 

C  PM(K)  MEAN  PRESSURES 

0  PV(K)  STANDARD  DEVIATION  OF  P(T) 

C  PCER(K)        60-SECOND  LOADING  OF  EQN.  20 

e  PDAG(K)        60-SECOND  LOADING  OF  EQN.  21 

C  PF(K)  #  OF  FAILURES  FOR  P(T)/#  OF  PANELS  TESTED 

C         UNITS  OF  PRESSURE  VALUES  =  PSF 

C  VELOCITY  IN  FT/ SEC 

C 

C 

DO  5556KK=1 ,NSTORM 
5556    WRITE (6, 5) UMEAN ( KK ) , PM( KK ) , PV ( KK ) , PCER ( KK ) , PDAG ( KK ) , PF ( KK ) 
GO  TO  5559 
5558  WRITE(6,9) 

DO  55601= 1,NST0RM 
P60(I)  =  1HJ».»P60(I) 

5560  WRITE(6,112)P60(I),PCER(I),PDAG(I),PF(I) 

5561  CONTINUE 
C 

6000  FORMAT(1X,5F10.4) 

9000    FORMATdX,'  THE  VALUE  OF  ILOAD,  ISTREN,  OR  ISTRES  IS  WRONG') 

2  F0RMAT(1X,I5,6X,I5,7X,I5,5X,I5,2X,F12.4) 

3  FORMATCIX, '  FAILURE  #  '.IX,'  POSITION  ',1X,'  DIRECTION  ', 
•IX,'  TIME     '.IX.'S  INITIAL  ') 

5  FORMAT(1X,6F10.5) 
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6  F0RMAT(3X,'  MEAN  U  ',3X,'  MEAN  P  ',1X,'  STAN  P'.IX, 

PI       ',1X,'  P2     ' ,1X, 'PF(PMEAN)  ') 

7  FORMATdX,'  STORM  NUMBER  ',110) 

8  F0RMAT(16C) 

9  F0RMAT(3X, '  PRESSURE  ',5X,'  PI    '.IX, •  P2    •  ,7X, 
»'  PF  ') 

112  FORMAT(1X,F10.2,5X,F10.2,5X,F10.5,5X,F10.5) 

kH        FORMATdX, 6F10. 4) 

C 

C 

5559  CONTINUE 
STOP 
END 

C 

SUBROUTINE  DIR(RSIG) 
C  FIND  DIRECTIONAL  MULTIPLIERS  FOR  STRESSES  AT  VARIOUS 

C  ANGLES... 
C 
C 

DOUBLE  PRECISION  PH,RSIG(36 ,6) ,C(6) ,SII(6) ,S2S1 (36) 
DATA  S2S1/ 1.0,. 98,. 94,. 89,. 79,. 68,. 52,. 12,. 97,. 93, 
*. 86,. 77,. 65,. 47,. 05,. 89,. 82,. 71,. 57,. 36, -.09,. 73,. 61, 
•.45,. 21, -.24,. 47,. 28,. 02, -.38,. 83, -.17, -.52, -.39, -.66, 
»-.8l/ 
DO  5K=1,6 

PH=(K-1 )*3. 14159/10. 

C{K)=DCOS(Ph)«DCOS(PH) 
5  SII(K)=DSIN(PH)»DSIN(PH) 

DO  10J=1,36 
DO  11K=1,6 
11  RSIG(J,K)=C(K)  +  S2S1(J)»SII(K) 

10  CONTINUE 
RETURN 
END 

C 

SUBROUTINE  PEQ(NDUR,DT,RN,P,P1 ,P2) 

C 
C 

C  THIS  ROUTINE  CALCULATES  THE  60-  SECOND  LOADINGS  OF 

C  EQNS.  20  AND  21 

C 

C 

C 

DOUBLE  PRECISION  RN ,R00T2, PEFF( 2) , PMAX, P(») 
•,AREA,P1,P2,PD(900) 
R00T2=1./RN 
PMAX=0.0 
IlrNDUR 
DO  10001=1,11 
PMAX=DMAX1(PMAX,P(I)) 
1000  CONTINUE 

DO  2000K=1,2 

PEFF( 1 ) = ( ( DT/60 . ) *»R00T2 ) *PMAX 

IF{K  .EQ.   1)  GO  TO  3000 

DO  25001= 1,NDUR 
2500  PD(I)=P(I)»»RN 
C 
C 

C  INT  IS  AN  INTEGRATION  ROUTINE 

C 

C 

CALL  INT ( PD , NDUR , 1 . 0 , AREA ) 

C 
C 

PEFF( 2) = ( AREA/60 . ) »»R00T2 
P2=PEFF(2) 
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3000  P1=PEFF(1) 
2000  CONTINUE 

RETURN 

END 

C 
C 

SUBROUTINE  INT(P,NDUR,H,AREA) 

C 
C 

C  TRAPEZOIDAL  RULE 
C 

DOUBLE  PRECISION    AREA, P( 900) 
N=NDUR-1 
SUM=0.0 
DO  501=2, N 
50  SUM=SUM+P(I) 

AREA=0.5»H»(  P(1)  +  2.0»SUM  +  P(NDUR)  ) 

RETURN 

END 

C 
C 

SUBROUTINE  PWRITE ( IT , MRI , DT , CP , PMEAN , P , NDUR ) 

C 
C 

C  SUBROUTINE  TO  PRINT  INPUT  PRESSURE  SIMULATION  PARAMETERS 

C 

C 

DOUBLE  PRECISION  P(900) 

REAL  IT 

WFITE(6,1) 

WRITE(6,2)IT,MRI,DT,CP 
WRITE (6, 3) NDUR 
WRITE (6, 30) PMEAN 
WRITE(6,4)(  P(I),  1=1, NDUR  ) 
WRITE(6,1) 

1  FORMAT(IHI) 

2  F0RMAT(//,1X,'  LOADING  SIMULATION  ',/, 
•5X,'  INPUT  MEAN  VELOCITY  (FT/SEC)  =',F10,4,/, 
»5X, •  MEAN  RECURRENCE  INTERVAL  (YEARS)  =',110,/, 
»5X, '  TIME  STEP  IN  SECONDS  =',F10.4,/, 
*5X, '  PRESSURE  COEFFICIENT,  CP  =',F10.4) 

3  F0RMAT(5X,'  P(T),  T=1,   ',15,'  SECONDS   ') 

30        F0RMAT(6X,'  MEAN  PRESSURE  IN  PSI  =  ',F10.4) 

4  FORMAT (5F1 0.4) 
RETURN 

END 

C 

SUBROUTINE  GWRITE ( SIDEA , TH , E , PR , A , KIC , NPANEL , RN , BP ) 

C 
C 

C  SUBROUTINE  TO  PRINT  INPUT  PLATE  PARAMETERS... 
C 

c 

DOUBLE  PRECISION    KIC, A, PR, E,RN, BP, SIDEA 

WRITE(6,1) 

WRITE (6, 2) NPANEL 

WRITE(6,3)SIDEA,TH,E,PR,A,KIC,RN,BP 

1  FORMAT(IHI) 

2  FORMATdX,'  ••  GLASS  PANEL  CHARACTERISTICS  • 
•2X,I5,'  PANELS  ARE  TO-BE  TESTED  ',/) 


FORMAT (/,5X, 

•  5X, 

•  5X, 

•  5X, 

•  5X, 

•  5X, 


LENGTH  OF  EACH  SIDE 
THICKNESS 

MODULUS  OF  ELASTICITY 
POISSONS  RATIO 
VALUE  OF  A 

STRESS  INTENSITY  FACTOR 


II, 

,F10.4,/, 
,F10.4,/, 
,E10.4,/, 
,F10.i»,/, 
,F10.4,/, 
,F10.4,/, 
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•  5X,'  POWER  N  IN  EQUATION  =',F10.1»,/, 

•  5X, '   1/B    IS  CALCULATED  TO  BE  = ' , F1 0.4 , /// ) 

RETURN 
END 

C 

SUBROUTINE  TWRITE ( ILOAD , ISTRES , ISTREN , NDUR ) 

C 

c 
c 

C  PROGRAM  TO  WRITE  TEST  PARAMETERS 

C 

C 

C 

WRITE(6,1) 

IF  (  ILOAD  .EQ.   1  )WRITE(6,2) 
IF  (  ILOAD  .EQ.  2  )WRITE(6,1») 
IF  (  ILOAD  .GT.  3  )WRITE(6,31) 
IF  (  ILOAD  .EQ.  3  )WRITE(6,32) 
WRITE (6, 3) NDUR 
WRITE(6,5) 

IF  (  ISTREN  .EQ.   1  )WRITE(6,6) 
IF(  ISTREN  .EQ.  2  )WRITE(6,7) 
WRITE(6,8) 

IF(  ISTRES  .EQ.  1  )WRITE(6,9) 
IF  (  ISTRES  .EQ.  2  )WRITE(6,10) 
WRITE (6, 3) NDUR 


1 

FORMATdX, ' 

TEST  PARAMETERS  *•»',//, 

•2X, '  #1 

PRESSURE  LOADING  ',/) 

2 

F0RMAT(7X, ' 

SIMULATED  P(T)  ') 

3 

FORMAT (7X, • 

FOR  ' ,15, '  SECONDS  ' ,/) 
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F0RMAT(7X, • 

NOT  CONSIDERED  ') 

32 

F0RMAT(7X, • 

READ  IN  FROM  FILE  ') 

4 

F0RMAT(7X, • 

CONSTANT  PRESSURE  ') 

5 

F0RMAT(2X, • 

#2          INITIAL  STRENGTH  DESCRIPTION' 

6 

F0RMAT(7X, ' 

TO  BE  GENERATED  FROM  A  WEIBULL  DISTRIBUTION  •) 

7 

F0RMAT(7X, • 

CONSTANT  VALUE  ') 

8 

FORMAT (2X, ' 

#3  ','  STRESS  TIME  HISTORY  ',/) 

9 

F0RMAT(7X, ' 

GENERATED  FROM  THE  PRESSURE  LOADING 

') 

10 

F0RMAT(7X, ' 

CONSTANT  SIGMA  VALUE  •) 

RETURN 
END 

C 

SUBROU TINE  PRES ( R 1 , R2 , NDI V , DS , DT , NDUR , Z , ZO , IT , CP , PM , PVAR , P ) 

C 

c 

C  THIS  SUBROUTINE  COMPUTES  THE  QUASI-STATIC  LOADING 

C  AT  THE  STAGNATION  POINT  USING  THE  ARIMA-COMBINED 

C  TECHNIQUE  AND  THE    SPECTRAL  EXPRESSION  IN  EQN.  18 

C 

C 

C 

C 

C 

c 

REAL  AO.IT.DENSTY.PJ ,Z,SUM,K1 ,C,UF 
DOUBLE  PRECISION  P(900),DS 
REAL  U(900),NORM(900) 

REAL  UFS , DT , PH1 1 , PHI2 , ACOV ( 3 ) , ACF( 3 ) , VARU ,  VARAU 
PI=3. 14159 
DENSTY=. 00258 

XMAX=R2 

XMINrRI 

POWER=5./3. 

DX= ( XMAX-XMIN ) /NDIV 

C 

c 


30 


c 

C  DEFINITIONS 


C 
(: 

C  IT=  MEAN  VELOCITY 

C  Z=ELEVAT1()N 

C  2.0=FRICTION  PARAMETER 

C  UF=FRICTION  VELOCITY 

C  ACOVrAUTOCOVARIANCE,  ACF=AUTOCORRELATION 

C  VARU,VARAU=  SCALING  VARIANCES  FOR  VELOCITY 

C  U=GENERATED  VELOCITY,  P=GENERATED  PRESSURES 
C 
C 


C  GGNML  IS  A  PSEUDO-RANDOM  #  GENERATOR  FROM  IMSL 

C 

C 

C  FOR  AN  EXPLANATION  OF  THE  METHOD,  SEE  "AUTOREGRESSIVE 

C  MODELING  OF  LONGITUDINAL  TURBULENCE  SPECTRA",  J.  OF 

C  INDUSTRIAL  AERODYNAMICS  AND  WIND  ENGINEERING,  1982, 

C  BY  REED  AND  SCANLAN. 

C 

C 

c 

A0r50.»Z,/IT 

UF=(0.4»IT)/(AL0G(12.»Z/Z0)) 

UFS=UF*UF 

C=(200.»UFS»Z)/IT 

C 

DO  10J=1,3 
1=1 
X=0.0 
SUM=1.0 
TAU=(2.«J)-2. 
TAU=PI»TAU 
2  X=X+DX 

FREQ=TAU»X 

FX=C0S(FREQ)/((1 .+A0»X)«»POWER) 

G=2.»FX 

SUM=SUM4-G 

1=1+1 

IF(I  .LT.  NDIV)GO  TO  2 

FB=COS ( FREQ ) / ( ( 1 . +AO»XMAX) »»POWER ) 

SUM=SUM+FB 

AREA=SUM«(DX/2.) 
10  AC0V(J)=(C/2.)*AREA 

DO  20J=1,2 
20  ACF(J)=AC0V(J+1)/AC0V(1) 
C 

D=1.-(ACF(1)»»2.) 
PHI1r(ACF(1)*(1 .-ACF(2)))/D 
PHI2=(ACF(2)-(ACF(1)**2.))/D 

WRITE(6,»)'  AUTOREGRESSIVE  MODEL  PARAMETERS  ARE  ',PHI1,PHI2 
VARU=AC0V(1) 

VARAU=(((1.+PHI2)«((1 .-PHI2)«»2.-(PHI1«»2. )))/(1 .-PHI2)) 
1  "VARU 

SDEVA=SQRT(VARAU) 
DO  30J=1,2 
30  U(J)=0. 

C 
C 
C 

CALL  GGNML(DS,NDUR,NORM) 

C 

c 

C  SET  UP  GENERATION  EQUATION  FOR  FLUCTUATING  VELOCITY 


31 


C  DEFINE  PMEAN 

PMEANr { 1 . / 1 41 . ) . 5»IT*IT»CP»DENSTY 

RCON=PMEAN/IT 

C 

DO  50J=3,NDUR 

50  U(J)=PHI1»U(J-1)  +  PHI2«U(J-2)  +  SDEVA«NORM( J) 

DO  60J=1,NDUR 
60  P(J)=PMEAN  +  RCON»U(J) 

C 
C 
C 

C  THIS  P  REPRESENTS  THE  TOTAL  LOAD,  MEAN  +  FLUCTUATING  COMPONENT 

C  UNITS  OF  P  =  PSI 

C 

C 

C 

SUM=0.0 

DO  5J=1,NDUR 

5  SUM=SUM+P{J) 
PM=SUM/NDUR 
SUM=0.0 

DO  6J=1,NDUR 

6  SUM=SUM  +  ((  P(J)  -  PM  )«»2  ) 
PVAR=SQRT(  SUM/(NDUR-1)  ) 

C 

C  WRITE{8,70)(P(I),I=1,NDUR) 

C70  FORMAT  (5F10.i») 

C 

RETURN 
END 

C 

SUBROUTINE  STRESS( FLEX, SIDEA , TH , NDUR , P , SIGMA ) 

C 
C 

C  THIS  PROGRAM  CALCULATES  THE  STRESS  FACTORS  FOR  A  SQUARE 

C  PLATE  ,  GIVEN  A  PRESSURE  TIME  HISTORY 

C 

C 

DOUBLE  PRECISION    TH, FLEX, SIDE A, 

»R1(36),R2(36),R3(36),Ri»(36),S1(36),S2(36),S3(36),Si»(36) 
DOUBLE  PRECISION      P(900) , SIGMA(36 , 900) 
COEF= ( . 00689*FLEX) / ( SIDEA •SIDEA*TH ) 

C 
C 

C  INPUT  INTERCEPT  AND  SLOPE  PARAMETERS  FOR  STRESS-PRESSURE 

C  LINEAR  RELATIONSHIPS 

C 

c 

DATA  R1/1H9.,151.,157.,163.,165.,157.,130.,78. J50.,155., 
•162., 165., 157., 131., 79., 151., 157., 162., 156., 13?., 82., 127., 
»1 27., 122., ion., 70,, 180., 181., 165., 120., 1 51., I4lt., 144., 

•144. ,143. ,157./ 
DATA  SI / . 058 , . 060 , . 066 , . 077 , . 092 , . 1 00 , . 094 , . 06 1 , . 059 , . 066 , 
*. 078,. 092,. 103,. 096,. 063,. 065,. 080,. 096,. 109,. 103,. 066, 
«. 086 , . 1 06 , . 1 1 9 , . 1 1 3 , . 07 1 , . 070 , . 084 , . 085 , . 056 , . 089 , . 086 , 
••063, .086, .080, .118/ 

C 

DATA  R2/  207. ,211. ,223. ,240. ,257. ,257. ,224. ,139. ,209., 
•221 . ,239. ,258. ,258. ,227. , 142. ,216. ,237. ,258. ,265. , 
•235., 148., 21 3., 233., 241., 217., 141., 250., 265., 250., 176., 
•238. ,230., 186. ,237. ,223. ,246./ 

DATA  S2/ . 034 , . 036 , . 036 , . 043 , . 054 , . 067 , . 067 , . 047 , 
•.037,. 039,. 045,. 055,. 068,. 069,. 048,. 044,. 048,. 059, 
•. 072 , . 075 , . 053 , . 059 , . 07 1 , . 089 , . 095 , . 068 , . 062, . 087 , 
••099,. 073,. 100,. 118,. 090,. 123,. 105,. 114/ 

C 
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DATA  R3/310.  ,319.  ,333.  ,370.  ,420.  ,i<57.  ,426.  ,279.  ,320. , 
•337.,  374., H25., 463.,  435.,  286.,  349.,  381., 436., 1)82., 
•460. ,307. , 392. ,445. ,507. ,502. ,345. ,436. ,525. ,549. , 
•396. ,539. ,583. ,457. ,607. ,607. ,587./ 

DATA  S3/ . 0 1 9 , . 0 1 8 , . 0 1 9 , . 023 , . 032 , . 045 , . 053 , . 038 , . 022 , 
•.023,. 025,. 034,. 047,. 053,. 039,. 028,. 031,. 037,. 049,. 055, 
•.041 ,.038,. 042,. 052,. 060,. 046,. 050,. 057,. 070,. 058,. 058, 
•.083,. 077,. 096,. 082,. 105/ 

DATA  R4/407., 410., 430., 484., 581., 687., 689., 471., 428. , 
•454., 500., 593., 697., 700., 480. ,490. ,538. ,620. ,725. ,735. , 
•51 1. ,584., 656. ,769. ,803. ,576. ,688. ,812. ,901 . ,688. ,827. , 
•998. ,843., 1086. , 1016. , 1114./ 

DATA  S4/ . 029 , . 029 , . 027 , . 026 , . 032 , . 042 , . 049 , . 036 , . 030 , . 03 1 , 
». 029,. 032,. 041,. 048,. 037,. 033,. 034, .034,. 042,. 049,. 037, 
••037,. 038,. 043,. 050,. 039,. 044,. 046,. 056,. 048,. 038,. 064, 
•.064, .072, .085,. 110/ 


DO  200I=1,NDUR 
IF(P(1).LE.1000)GO  TO  3 

IF(P(I)   .GT.   1000  .AND.  P(I)   .LE.  2000)  GO  TO  4 
IF(  P(I)     .GT.  2000  .AND.   P(I)   .IF.   5000)  GO  TO  5 
IF(P(I)     .GT,   5000  .AND.   P(I)   .LE.   10000)  GO  TO  6 
IF(P(I)     .GT.  10000)  GO  TO  7 


3  CONTINUE 

SIGMA(1,I)=  COEF«(202  -  DLOG(P(I))»(  99.6  -  13. 3*DL0G(P(I) ) ) ) 

SIGMA(9,I)=SIGMA(1 ,1) 

SIGMA(16,I)=SIGMA(1,I) 

SIGMA (22, I)=SIGMA( 1,1) 

SIGMA(2,I)=SIGMA(1 ,1) 

SIGMA(3,I)=SIGMA(2,I) 

SIGMA(10,I)=SIGMA(2,I) 

SIGMA(4,I)=C0EF»(267  -  DLOG(P(I) )«( I3I  -16.9*DL0G(P(I) ) )  ) 

SIGMAdI  ,I)  =  SIGMA(4,I) 

SIGMA(17,I)=SIGMA(4,I) 

SIGMA(5,I)  =  C0EF»(441  -  DL0G(P(I))^(  201  -  23- 5^DL0G( P(I) ) ) ) 

SIGMA(12,I)=SIGMA(5,I) 

SIGMA(18,I)=SIGMA(5,I) 

SIGMA(23,I)=SIGMA(5,I) 

SIGMA(6,I)=  C0EF»(  676-  DLOG(P(I) )^(292-31 .6*DL0G{P(I) ) ) ) 

SIGMA(13,I)=SIGMA(6,I) 

SIGMA(19,I)=SIGMA(6,I) 

SIGMA(24,I)=SIGMA(6,I) 

SIGMA(28,I)=SIGMA(6,I) 

SIGMA(7,I)=C0EF»(  4.96  +  .130^P(I)  ) 

SIGMA(14,I)=SIGMA(7,I) 

SIGMA(20,I)=SIGMA(7,I) 

SIGMA(25,I)=  C0EF»(27.4  +  .0790*P(I)  ) 

SIGMA(29,I)=C0EF»(  30.2  +  .0885*P(I)  ) 

SIGMA(32,I)=  C0EF»(  36.1  +  .0969^P(I)  ) 

SIGMA(8,I)=  COEF»(  1.93  +  .0785^P(I)) 

SIGMA(15,I)=SIGMA(8,I) 

SIGMA(21,I)=SIGMA(8,I) 

SIGMA(26,I)=  COEF»(  2.70  +  .0384»P(I)) 

SIGMA(30,I)=  C0EF*(  6.52  +  0.109^P(I)) 

SIGMA(33,I)=  C0EF»(  8.37  +  .120»P(I)) 

SIGMA(35,I)=  COEF*(  7.93  +  .142^P(I)) 

SIGMA(36,I)=  C0EF»(  7.07  +  0.156^P(I)) 

SIGMA ( 27 , I ) =SIGMA ( 36 , I ) 

SIGMA(31,I)=SIGMA(36,I) 

SIGMA ( 34 , I ) = SIGMA ( 36 , I ) 
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c 
c 

GO  TO  200 

C 

c 

k  CONTINUE 
DO  42J=1,36 

42        SIGMA(J,I)=  COEF»(  R1(J)  +  S1(J)»P(I)  ) 

C 

C 

GO  TO  200 

C 
C 

5  CONTINUE 

DO  52J=1,36 

52         SIGMA(J,I)=COEF*(  R2(J)  +  S2(J)«P(I)  ) 
GO  TO  200 

C 

6  CONTINUE 

DO  62J=1,36 

62  SIGMA(J,I)=  C0EF«(R3(J)  +  S3(J)»P(I)  ) 
C 

c 

GO  TO  200 

C 

c 

7  CONTINUE 

DO  72J=1,36 

72  SIGMA(J,I)=  C0EF«(R4(J)  +  S1»(J)»P(I)  ) 
C 

200  CONTINUE 
C 

RETURN 
END 

$BEND 
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Table  1.     Estimated  Probabilities  of  Failure  Based  on  Stress  Time 
History,  and  on  l-Minute  Loads  Estimated  in  Accordance 
with  References  12  and  3.     Units  of  pressure  =  psf. 


p 

C 
P60 

D 
P60 

Pf  (P) 

Pf(P60  -  P60^ 

Pf(P60  -  P60) 

10 

10.8 

12.6 

0.001 

0.002 

0.003 

15 

16.2 

18.9 

0.015 

0.004 

0.008 

20 

21.6 

25.2 

0.050 

0.012 

0.022 

25 

27.0 

31.5 

0. 150 

0.035 

0.069 

30 

32.4 

37.8 

0.400 

0.084 

0.160 
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Principal 

stress 

directions 


Figure  3.     Schematic  representation  of  angles  a]^  =  kAa 
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Figure  4.     Cumulative  distribution  function  for  p(t) 
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0.20 


TME,  t  (sj 

Figure  5.     Evolution  of  strength  and  stresses  with  time  for  breaking  panel 
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